##########################################################################
## To create simulation of a multinomial choice with social interactions##
## We take true (N,nei) distribution from our data, assuming those belonging to ##
## the same street interact with each other ##
## Created  07/08/2014
## Modified 16/07/2019
## Where we set L=5 as in our original submission
## Where we guarantee that KS probabilities are within [0,1]
## We guarantee inner loop (fixed point expectations) convergence using tol 10e-8 for the abs-max stopping criterion
## We follow KS 2012 and KS 2018 using outer loope iterations k=50 and guaranteeing inner loop convergence (see above)
## We set coefficients to be equal to the original submission
## by: Jose Alberto Guerra
##########################################################################

rm(list=ls(all=TRUE))

#########################
####    Preamble     ####
#########################

# Definig Globals 
#Root <- "/Users/josealbertoguerra"
Root <- "C:/Users/ja.guerra/Dropbox"
WD <- paste(Root,"/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes",sep="")
Codes <- paste(WD,"/Submission/Final_ReStat/Codes/C_ANALYSIS",sep="")
#Codes <- paste(WD,"/Codes/RCodes/EEACodes",sep="")
Results <- paste(WD, "/Results/ReStat/",sep="")
#Results <- paste(WD, "/Results/EEA/",sep="")
setwd(WD)
Graphs<- paste(getwd(),"/Results/ReStat/Graphs/",sep="")


source(paste(Codes,"/Mlogit_MontecarloFunctions.R",sep=""))
#  Then try/install packages...
foo(c("spdep","RANN","gumbel","fExtremes","mlogit","mnlogit","nnet","reshape2","grDevices","data.table","ggplot2","MASS","dplyr")) 


# Number of members per parish $ n_g \in \left\{100 300 500 right\} $ 
# should be changed below and xx in line "Ngs <- c(xx)"


.pardefault<-par() # defining default parameters graphs


#########################
### Global definitions###
#########################
# Seed for the parameters
set.seed(12345)
#Load data with Global Coefficients
load(paste(Codes,"/datatrueJGood.RData",sep="")) #datatrueJ
load(paste(Codes,"/datatrueaGood.RData",sep="")) #datatruea
load(paste(Codes,"/datatruedGood.RData",sep="")) #datatrued
load(paste(Codes,"/datatruea0Good.RData",sep="")) #datatruea0

# Parameters
L <- 5 # number of choices
K <- 1 #n indi covariates
S <- K #n group covariates
tresh <- 4*(1-1/L) #threshold I found for unique equilibria if all coeff are equal
tresh <- tresh+.5
Jtresh <- 0 #with 1/4 probability J<0
Jcheck <- if (L>2) (runif(L,0,1)<Jtresh) else (runif(L-1,0,1)<Jtresh)

#The number of worlds which will give us the distribution on the estimated parameters
Nworlds <- 100

#The number of parishes per world
#Nparishes <- c(5,10,20,50) #if we iterate 
Nparishes <-  c(150)# $ \# G$c(100,150)
Ngs <- c(500) #Number of obs per group $ n_g $ c(100,300,500)

for(Ng in Ngs) {
  for(Nparish in Nparishes) { 
    sink(paste(getwd(),"/Results/ReStat/KS/Montecarlo",Ng,"ViPMLNPL",Nparish,".txt",sep="")) ###
    print("##############")
    print(paste("Ng: ",Ng,sep=""))
    print("##############")
    print(paste("World with N parishes: ",Nparish,sep=""))
    
    #Outer interation, to keep track of the parishes
     maxitt0 <- 2000 #maximum number of iterations
    # 
     a0ea0 <- as.matrix(matrix(rep(0,maxitt0*(L-1)*4),nrow=maxitt0))
     a0ea0ks <- as.matrix(matrix(rep(0,maxitt0*(L-1)*4),nrow=maxitt0))
    # truea0 <- as.matrix(matrix(rep(0,maxitt0*L),nrow=maxitt0))
     aea <- as.matrix(matrix(rep(0,maxitt0*(L-1)*4),nrow=maxitt0))
     aeaks <- as.matrix(matrix(rep(0,maxitt0*(L-1)*4),nrow=maxitt0))
    # truea <- truea0
     ded <- as.matrix(matrix(rep(0,maxitt0*(L-1)*4),nrow=maxitt0))
     dedks <- as.matrix(matrix(rep(0,maxitt0*(L-1)*4),nrow=maxitt0))
    # trued <- truea0
     JeJ <- as.matrix(matrix(rep(0,maxitt0*(L)*4),nrow=maxitt0))
     JeJks <- as.matrix(matrix(rep(0,maxitt0*(L)*4),nrow=maxitt0))
    # trueJ <- truea0
    nn0 <- as.matrix(matrix(rep(0,maxitt0),nrow=maxitt0))
    nnei0 <- as.matrix(matrix(rep(0,maxitt0),nrow=maxitt0))
    nn0ks <- as.matrix(matrix(rep(0,maxitt0),nrow=maxitt0))
    nnei0ks <- as.matrix(matrix(rep(0,maxitt0),nrow=maxitt0))
    # datalistuue <- vector(length=maxitt0,mode="list") #keep track of correlated effect
    # 
    choices <- as.matrix(matrix(rep(0,maxitt0*(L)),nrow=maxitt0))
    choicesks <- as.matrix(matrix(rep(0,maxitt0*(L)),nrow=maxitt0))
    #Define coefficients, we will use those from 5th iteration
    #a0
    a0 <- as.matrix(matrix(unlist(datatruea0[5,c(1:L)]),nrow=1))
    #a
    a <- as.matrix(matrix(unlist(datatruea[5,c(1:L)]),nrow=1))
    #d
    d <- as.matrix(matrix(unlist(datatrued[5,c(1:L)]),nrow=1))
    #J
    J <- as.matrix(matrix(unlist(datatrueJ[5,c(1:L)]),nrow=1))
    #Correlated effect
    sig <- c(0.13 ,0.08, 0.18, 0.05 ,0.1)
    useed <- sapply(rep(1,L),function(x) rep(rnorm(x),1))*sig
    
    
    nogooditt0 <- 0 #keep track of those worlds for which we don't have enough variation in choices 
    gooditt0 <- 0 
    itt0 <-1
    while (itt0 <=maxitt0 & gooditt0 < Nworlds ) {
      print(paste("**** Outer worlds iteration: ",itt0,sep=""))
      #Simulations for maxitt parishes
      maxitt <-400
      coordslist <- vector(length=maxitt,mode="list") #keep track of coordinates per group
      truewlist <- vector(length=maxitt,mode="list") #keep track of the true spatial weights per group
      wklist <- vector(length=maxitt,mode="list") #keep track of another set of spatial weights per group
      datalist <- vector(length=maxitt,mode="list") #keep track of the data
      nn <- as.matrix(matrix(rep(0,maxitt),nrow=maxitt)) #keep track of the numer of observations per group
      nnei <- as.matrix(matrix(rep(0,maxitt),nrow=maxitt)) #keep track of the number of neighbours per group
      uu <- as.matrix(matrix(rep(0,maxitt*L),nrow=maxitt)) #keep track of the group fixed effect
      
      nogooditt <- 0 #keep track of those parishes for which we don't have all choices 
      gooditt <- 0 
      itt <- 1
      while (itt <= maxitt & gooditt < Nparish ) {
        print(paste("  Inner iteration: ",itt,sep=""))
        print(paste("  Spatial dimension: ",itt,sep=""))
        #########################
        ### Spatial dimension ###
        #########################
        # Number of parish members and min neigh
        #minneiN <- NAPPsocials[round(runif(1,1,nrow(NAPPsocials)),0),c("minnei", "N")] # NAPPsocials[round(runif(1,1,nrow(NAPPsocials)),0),c("mnei", "N")] #
        #n <- minneiN[[2]]
        #nei <- minneiN[[1]]
        n <- Ng #Number of obs per group
        nei <- 10 #number of minimum nei
        nn[itt,1] <- n
        nnei[itt,1] <- nei
        # Coordinates
        coords <- cbind(runif(n,0,runif(1,1,2)), runif(n,0,runif(1,1,2)))
        coordslist[[itt]] <- coords
        # Distance Matrix
        #D <- as.matrix(dist(coords))
        kn<-10 #ceiling(n*.05)
        
        # Nearest neighbour
        kneighbours1<- knn2nb(knearneigh(coords,nei))
        #k-nearest neighbour
        kneighbours<-knn2nb(knearneigh(coords,kn))
        wk<-nb2listw(kneighbours, glist=NULL, style="W", zero.policy=NULL) #R row standardised
        wklist[[itt]] <- wk
        
        # Finding the smallest distance so each point is linked according to minnei
        all.linked <- max(unlist(nbdists(kneighbours1, coords)))
        dneighbours<-dnearneigh(coords,0,all.linked)
        wd<-nb2listw(dneighbours, glist=NULL, style="W", zero.policy=NULL)
        rm(kneighbours1)
        
        #True neighbours
        truew <- wd #or wk
        truewlist[[itt]] <- truew
        # Graphs Maps
        #plot(coords,col="blue")
        #points(coords[7,1], coords[7,2], pch=18)
        #par(mfcol=c(1,2), pty="s")
        #plot(kneighbours,coords, col="blue")
        #if (itt %in% c(4,27,34,43,50,62,78,76,90)) {
        #cairo_ps(paste("SampleMapNeigh_n",n,".eps",sep=""), width=7, height=7) #to keep transparency from the intervals
        #print(plot(dneighbours,coords, col="red"))
        #title(paste("n obs:",n,". av n peer:",round(mean(unlist(lapply(dneighbours,length))),2)))
        #dev.off()
        #}
        #par(.pardefault)
        
        #########################
        #### Data dimension  ####
        #########################
        
        ###Group specifics
        print(paste("  Data dimension ",itt,sep=""))
        # Covariates
        c <- as.matrix(rep(1,n))
        #X <- as.matrix(c(rnorm(n,1,1.5)))
        X <- as.matrix(c(rnorm(n,runif(1,.8,1.2),runif(1,.2,2.2)))) 
        Xw <- as.matrix(lag.listw(truew,X))
        #Y <- as.matrix(matrix(runif(n*L,0,1),nrow=n)) #Y <- as.matrix(matrix(rnorm(n*L),nrow=n)) # to test whether non-identification is due to linear combinaiton or lack of sufficient support in expec
        #Y <- Y/apply(Y,1,sum)
        
        #correlated effect
        u <- t(apply(as.matrix(sapply(rep(1,L),function(x) rep(rnorm(x),n))),1,function(x) x*sig))
        #u <- as.matrix(matrix(rep(useed,n),ncol=L,byrow=TRUE))
        uu[itt,] <- apply(u,2,mean)
        
        # pre-Utility and initial beliefs educated guess
        epsilon <- matrix(rgev(n*L, xi = 0, mu = 0, beta = 1),ncol=L)
        perchoice <- runif(L) #to give some starting value, however it seems not to matter given our fixed point iteration. To see this replace: expec <- as.matrix(matrix(rep(0,n*L),ncol=L))
        perchoice <- perchoice/sum(perchoice)
        expec <- t(rmultinom(n,1,perchoice)) 
        expec <- exp(c%*%a0 + X%*%a + Xw%*%d + t(apply(as.matrix(lag.listw(truew,expec)),1,function(x) x*J)) + u)/apply(exp(c%*%a0 + X%*%a + Xw%*%d +t(apply(as.matrix(lag.listw(truew,expec)),1,function(x) x*J)) + u),1,sum)
        
        # Iteration to find fixed point in rational heterogenous beliefs
        #definitions
        it <- 1 
        maxiter <- 300
        expectit <- matrix(rep(0,maxiter*L),nrow=maxiter)
        ndexpectit <- L
        zero <- 1E-8
        #iteration
        while (it<=maxiter & ndexpectit > zero) {
          #print(it)
          expectit[it,] <-apply(expec,2,mean)
          expec1 <- exp(c%*%a0 + X%*%a + Xw%*%d+t(apply(as.matrix(lag.listw(truew,expec)),1,function(x) x*J)) + u)/apply(exp(c%*%a0 + X%*%a + Xw%*%d+t(apply(as.matrix(lag.listw(truew,expec)),1,function(x) x*J)) + u),1,sum)
          if (it>=2) ndexpectit <- sum(distr(expec,expec1))
          expec <- expec1
          it <- it+1
        }
        if (it<=maxiter & ndexpectit < zero) print(paste("Fixed point in beliefs convergence reached at iter: ",it,sep="")) else {
          print(paste("Fixed point in beliefs FAILED to converge after it: ",it,sep=""))
          nogooditt <- c(nogooditt,itt) 
        }
        
        
        # Utility,weighted heterogeneous beliefs and choices
        #weighted Heterogenous beliefs
        expecw <- as.matrix(lag.listw(truew,expec))
        #expecwk <- as.matrix(lag.listw(wk,Y)) #to test whether non-identification is due to linear combination or lack of sufficient support in expec, it seems to be the case
        
        #utility
        V <- c%*%a0 + X%*%a + Xw%*%d + t(apply(expecw,1,function(x) x*J))+ epsilon + u 
        #V <- c%*%a0 + X%*%a  +  Xw%*%d + t(apply(Y,1,function(x) x*J))+epsilon
        
        #choices
        choice <- as.matrix(max.col(V))
        if (length(unique(choice)) < L)  nogooditt <- c(nogooditt,itt) 
        apply(expec,2,mean)
        
        #additional definitions to be exported in the data
        s <- as.matrix(apply(V,2,function(x) (x==as.matrix(apply(V,1,max)))*1))
        swk <- as.matrix(lag.listw(wk,s))
        swd <- as.matrix(lag.listw(wd,s))
        sw <-  as.matrix(lag.listw(truew,s))
        
        #Exporting data
        dataw <- data.frame(as.matrix(cbind(rep(itt,n),seq(1,n),X,Xw,sw,expecw,swk,choice)))
        colnames(dataw) <- c("PARID","ID","x","xw",paste("sw.",c(1:L),sep=""),paste("expw.",c(1:L),sep=""),paste("swk.",c(1:L),sep=""),"choice")
        datalist[[itt]] <- dataw
        
        itt <-itt+1
        nogooditt <- unique(nogooditt)
        gooditt <- itt - length(nogooditt)
        print(paste("good samples:",gooditt))
      }
      
      
      #Identifying those worlds for which we could not get enough parishes
      if (gooditt < Nparish ) nogooditt0 <- c(nogooditt0,itt0) else nogooditt0 <- nogooditt0
      
      #Doing the estimation for those worlds for which we have good variation
      if (!(itt0 %in% nogooditt0)) {
        
        #removing those groups that do not have variation in terms of choices
        if (length(nogooditt) >1) nogooditt <- c(nogooditt[2:length(nogooditt)],itt:maxitt) else nogooditt <- c(itt:maxitt)
        
        coordslist <- coordslist[-nogooditt]
        truewlist <- truewlist[-nogooditt]
        wklist <-  wklist[-nogooditt]
        datalist <- datalist[-nogooditt]
        nn <- nn[-nogooditt]
        nnei <- nnei[-nogooditt]
        uu <- uu[-nogooditt,]
        
        pooleddatalist <- rbindlist(datalist)
        print(table(pooleddatalist$choice)/sum(table(pooleddatalist$choice)))
        ltresh <- 0.05
        if (sum(table(pooleddatalist$choice)/sum(table(pooleddatalist$choice))>ltresh)==L) {
          #it there is enough variation in choices, plot
          if (!((itt-1) %in% nogooditt) & itt0 == 1 ) {
            #ploting network of the last parish 
            cairo_ps(paste(Graphs,"/SampleMapNeigh_Nparish",Nparish,"Ng",Ng,".eps",sep=""), width=7, height=7) #to keep transparency from the intervals
            #print(plot(dneighbours,coords, col="red", pch=c(choice)))
            print(plot(dneighbours,coords, col="red", pch=c(choice), lty=1))
            for (pc in c(1:L)) {
              if (length(coords[choice==pc,])/2>1) print(points(coords[choice==pc,],pch=pc,col="black")) else  print(points(coords[choice==pc,][1],coords[choice==pc,][2],pch=pc,col="black"))
            }
            
            #legend.loc <- c(bbox(coords)[1, 1]+.01 , bbox(coords)[2, 1]+.3 )
            legend("bottomleft", c("neigh", paste("c.",c(1:L),sep="")), col = c("red", rep("black",L)), text.col = "black", lty = c(1,rep(NA,L)), pch = c(NA, seq(1:L)),merge = TRUE, bg = "gray90",cex=.75)
            title(paste("n obs:",n,". av n peer:",round(mean(unlist(lapply(dneighbours,length))),2)))
            dev.off()
            
            #Plot convergence in beliefs
            colnames(expectit)<-c(paste("expec.",c(1:L),sep=""))
            dataexpectit<-as.data.frame(expectit[c(1:(it-1)),])
            dataexpectit$iter <- c(1:nrow(dataexpectit))
            meltdataexpectit<-melt(dataexpectit,id="iter")
            
            cairo_ps(paste(Graphs,"/SampleExpect_Nparish",Nparish,"Ng",Ng,".eps",sep=""), width=7, height=7) #to keep transparency from the intervals
            print(ggplot(meltdataexpectit,aes(x=iter,y=value,colour=variable,group=variable)) + geom_line()+ ylab(expression(paste(J[y]^{t})))+xlab("t"))
            dev.off()
          }
          
          #Poduce the long format and pooled data
          
          datalistlong <- lapply(datalist,function(x) tolong(x,"choice",c("sw.","expw.","swk.")))
          
          pooleddata <- rbindlist(datalistlong)
          
          pooleddataks <- rbindlist(datalistlong)
          
          #       save(pooleddata, file=paste("PooleddataCityGood",itt0,".RData",sep=""))
          #       save(datalist, file=paste("datalistGood",itt0,".RData",sep=""))
          #       save(coordslist, file=paste("coordslistGood",itt0,".RData",sep=""))
          #       save(truewlist,file= paste("truewlistGood",itt0,".RData",sep=""))
          #       save(wklist, file= paste("wklistGood",itt0,".RData",sep=""))
          #       save(uu, file=paste("uuGood",itt0,".RData",sep=""))
          
          #pooleddataa<-load(paste("PooleddataCity",itt0,".RData",sep=""))
          rm(datalistlong)
          
          ##########################
          ####    Estimation    ####
          #### True expectations####
          ##########################
          
          print("Estimating...")
          print("Real J: ")
          print(J)
          #keeping track of the coefficients
          #True
          a0ea0[itt0,c((3*L-2):(4*(L-1)))] <- (a0-a0[1]+uu[1,]-uu[1,1])[2:L] ##given we control for PARID factors ce_l = c_l-c_0+u_{1,l}-u_{1,0}
          aea[itt0,c((3*L-2):(4*(L-1)))] <- (a-a[1])[2:L]
          ded[itt0,c((3*L-2):(4*(L-1)))] <- (d-d[1])[2:L]
          JeJ[itt0,c((3*L+1):(4*(L)))] <- J[1:L]
          nn0[itt0,1] <- sum(nn)
          nnei0[itt0,1] <- mean(nnei)
          choices[itt0,] <- table(pooleddata$alt*pooleddata$choice)[2:(L+1)]/sum(table(pooleddata$alt*pooleddata$choice)[2:(L+1)])
          #For KS
          a0ea0ks[itt0,c((3*L-2):(4*(L-1)))] <- (a0-a0[1]+uu[1,]-uu[1,1])[2:L] ##given we control for PARID factors ce_l = c_l-c_0+u_{1,l}-u_{1,0}
          aeaks[itt0,c((3*L-2):(4*(L-1)))] <- (a-a[1])[2:L]
          dedks[itt0,c((3*L-2):(4*(L-1)))] <- (d-d[1])[2:L]
          JeJks[itt0,c((3*L+1):(4*(L)))] <- J[1:L]
          nn0ks[itt0,1] <- sum(nn)
          nnei0ks[itt0,1] <- mean(nnei)
          choicesks[itt0,] <- table(pooleddataks$alt*pooleddataks$choice)[2:(L+1)]/sum(table(pooleddataks$alt*pooleddataks$choice)[2:(L+1)])
          
          fformula <- mFormula(choice ~ 1 | x + xw + factor(PARID)| expw)
    
          ######################
          #With observed expw##
          ######################
          
          print("1. Observed: expw")
          ml.m.w <-mnlogit(fformula, pooleddata, "alt")
          ml.m.wks <-mnlogit(fformula, pooleddataks, "alt") 
          #Estimated coefficientes
          a0e <- c(summary(ml.m.w)$coeff[paste("(Intercept):",c(2:L),sep="")])
          ae <- c(summary(ml.m.w)$coeff[paste("x:",c(2:L),sep="")])
          de <- c(summary(ml.m.w)$coeff[paste("xw:",c(2:L),sep="")])
          Je <- summary(ml.m.w)$coeff[paste("expw:",c(1:L),sep="")]
          print(Je)
          a0eks <- c(summary(ml.m.wks)$coeff[paste("(Intercept):",c(2:L),sep="")])
          aeks <- c(summary(ml.m.wks)$coeff[paste("x:",c(2:L),sep="")])
          deks <- c(summary(ml.m.wks)$coeff[paste("xw:",c(2:L),sep="")])
          Jeks <- summary(ml.m.wks)$coeff[paste("expw:",c(1:L),sep="")]
          print(Jeks)
          PARIDS <- c(unique(pooleddata$PARID))
          uenames<-paste(c(sapply(PARIDS, function(x) c(paste("factor(PARID)",x,sep="")), USE.NAMES=FALSE)),":",c(2:L),sep="")
          ue <- matrix(summary(ml.m.w)$coeff[c(sapply(sapply(PARIDS, function(x) c(paste("factor(PARID)",x,":",sep="")), USE.NAMES=FALSE),function(x) c(paste(x,c(2:L),sep="")), USE.NAMES=FALSE))],ncol=(L-1),byrow=TRUE)
          datauue <- as.data.frame(cbind(ue[2:gooditt,],t(apply((uu-uu[,1]),1,function(x) x-(uu-uu[,1])[1,]))[2:gooditt,2:L],nn[2:gooditt]))
          #notice ue_{p,l}=(u_{p,l} - u_{p,0})-(u_{1,l}-u_{1,0})
          colnames(datauue) <- c(paste("ue",c(2:L),sep=""),paste("u",c(2:L),sep=""),"nobs")
          
          #estimated
          a0ea0[itt0,c(1:(L-1))] <- a0e
          aea[itt0,c(1:(L-1))] <- ae
          ded[itt0,c(1:(L-1))] <- de
          JeJ[itt0,c(1:(L))] <- Je
          #For KS
          a0ea0ks[itt0,c(1:(L-1))] <- a0eks
          aeaks[itt0,c(1:(L-1))] <- aeks
          dedks[itt0,c(1:(L-1))] <- deks
          JeJks[itt0,c(1:(L))] <- Jeks
          
          ######################
          #Without observed expw##
          ######################
          print("2. Without Observed: expw")
          #wexercise <- c("w","wk") #w is true network, wk is just the first 10 neighbours
          wexercise <- c("w")
          datalist0 <- datalist
          pooleddatalist0 <- rbindlist(datalist0)
          pooleddata0 <- pooleddata
          coordslist0 <- coordslist
          truewlist0 <- truewlist
          wklist0 <- wklist
          uu0 <- uu
          #For KS
          datalist0ks <- datalist
          pooleddatalist0ks <- rbindlist(datalist0ks)
          pooleddata0ks <- pooleddataks
          coordslist0ks <- coordslist
          truewlist0ks <- truewlist
          wklist0ks <- wklist
          uu0ks <- uu
          #Iteration for each weighting matrix exercise
          for (wexer in wexercise) {
            
            pooleddata <- pooleddata0
            datalist <- datalist0
            pooleddatalist <- pooleddatalist0
            coordslist <- coordslist0
            truewlist <- truewlist0
            wlist <- truewlist0
            wklist <- wklist0
            uu <- uu0
            #For KS
            pooleddataks <- pooleddata0ks
            datalistks <- datalist0ks
            pooleddatalistks <- pooleddatalist0ks
            coordslistks <- coordslist0ks
            truewlistks <- truewlist0ks
            wlistks <- truewlist0ks
            wklistks <- wklist0ks
            uuks <- uu0ks
            weightslist <- list(wlist,wklist)
            
            #defining weight
            weigh <- weightslist[[which(wexercise==wexer)]]
            
            #General formula with intercept    
            fformula <-mFormula(as.formula(paste("choice ~ 1 | x + xw + factor(PARID)|", "s",wexer,sep="")))          
            print(paste("Observed: ","s",wexer,sep=""))
            
            #Outer Iteration: Maximum likelihood
            it0 <- 1 
            maxiter0 <- 50
            zero0 <- 1E-8
            Jeit <- matrix(rep(0,maxiter0*L),nrow=maxiter0) #Keep track of the coefficients
            a0eit <- matrix(rep(0,maxiter0*(L-1)),nrow=maxiter0) 
            aeit <- matrix(rep(0,maxiter0*(L-1)),nrow=maxiter0) 
            deit <- matrix(rep(0,maxiter0*(L-1)),nrow=maxiter0) 
            ndJeit <- L
            #For KS
            Jeitks <- matrix(rep(0,maxiter0*L),nrow=maxiter0) #Keep track of the coefficients
            a0eitks <- matrix(rep(0,maxiter0*(L-1)),nrow=maxiter0) 
            aeitks <- matrix(rep(0,maxiter0*(L-1)),nrow=maxiter0) 
            deitks <- matrix(rep(0,maxiter0*(L-1)),nrow=maxiter0) 
            ndJeitks <- L
            
            #To keep original pooleddata to implement Weidner (2012) method
            alter <- sort(unique(pooleddata$alt))
            pooleddata00 <- pooleddata0
            # To be able to use predict probabilities with mnlogit
            #pooleddata <- tolong(pooleddatalist0,"choice",c("sw.","expw.","swk."))
            pooleddata <- tolong(pooleddatalist0,"choice",c("sw.","expw.","swk."))
            pooleddataks <- tolong(pooleddatalist0ks,"choice",c("sw.","expw.","swk."))
            
            #to define outer iteration of expectations
            ndexpectito <- L
            ndexpectitoall <- matrix(0,nrow=maxiter0) #Keep track of the tolerance level
            expeco <- subset(pooleddatalist,select=paste("s",wexer,".",alter,sep=""))
            
            ndexpectitoks <- L
            ndexpectitoallks <- matrix(0,nrow=maxiter0) #Keep track of the tolerance level
            expecoks <- subset(pooleddatalistks,select=paste("s",wexer,".",alter,sep=""))
            
            while (it0<=maxiter0 & ndexpectito > zero0 & ndexpectitoks > zero0) {
              print(paste("ML iteration: ",it0, sep=""))
              #Estimating
              ml.m.w <-mnlogit(fformula, pooleddata, "alt")
              ml.m.wks <-mnlogit(fformula, pooleddataks, "alt")
              if(it0==1) { 
                ml.m.w0ks <- ml.m.wks
              }
              
              #Keeping coeff  
              Jeit[it0,] <- summary(ml.m.w)$coeff[paste("s",wexer,":",c(1:L),sep="")] 
              print(summary(ml.m.w)$coeff[paste("s",wexer,":",c(1:L),sep="")])
              deit[it0,] <- summary(ml.m.w)$coeff[paste("xw:",c(2:L),sep="")]
              aeit[it0,] <- summary(ml.m.w)$coeff[paste("x:",c(2:L),sep="")]
              a0eit[it0,] <- summary(ml.m.w)$coeff[paste("(Intercept):",c(2:L),sep="")]
              ndexpectitoall[it0,] <- ndexpectito 
              print(ndexpectitoall[it0,])
              #For KS
              Jeitks[it0,] <- summary(ml.m.wks)$coeff[paste("s",wexer,":",c(1:L),sep="")] 
              print(summary(ml.m.wks)$coeff[paste("s",wexer,":",c(1:L),sep="")])
              deitks[it0,] <- summary(ml.m.wks)$coeff[paste("xw:",c(2:L),sep="")]
              aeitks[it0,] <- summary(ml.m.wks)$coeff[paste("x:",c(2:L),sep="")]
              a0eitks[it0,] <- summary(ml.m.wks)$coeff[paste("(Intercept):",c(2:L),sep="")]
              ndexpectitoallks[it0,] <- ndexpectitoks 
              print(ndexpectitoallks[it0,])
              
              #Inner iteration, fixed point in beliefs
              it <- 1
              zero <- 1E-8
              #maxiter <- 1 #if iterating a la Aguirregabiria and Mira (2007) or a la Weidner (2012)
              maxiter <- 100 #If iterating with the fixed point a la Lee Lin Liu (2014) or Kasahara and Shimotsu (2018)
              expectit <- matrix(rep(0,maxiter*L),nrow=maxiter)
              expectitks <- matrix(rep(0,maxiter*L),nrow=maxiter)
              ndexpectit <- L
              ndexpectitks <- L
    
              #subsetting sample to 
              expec <- subset(pooleddatalist,select=paste("s",wexer,".",c(1:L),sep=""))
              expecks <- subset(pooleddatalistks,select=paste("s",wexer,".",c(1:L),sep=""))
              
              #iteration of beliefs
              while (it<=maxiter & ndexpectit > zero & ndexpectitks > zero) {
                print(paste("Inner iteration, fixed point in beliefs ",it,sep=""))  
                #keepint track of expectations
                expectitks[it,] <- apply(expecks,2,mean)
                
                alpha <- .1 #Kasahara and Shimotsu (2012)
                #Predict probabilities on the pooled data list
                tempprob <- ppredict.mnlogit(ml.m.w,pooleddata,probability=TRUE)
                tempprobks <- ppredict.mnlogit(ml.m.wks,pooleddataks,probability=TRUE)
                for(aa in sort(alter, decreasing=TRUE) ){
                  ap <- which(alter==aa)
                  
                  ##Aguirregabiria and Mira (2007)
                  if (it==1) {  
                    pooleddatalist[,paste("s",wexer,".",aa,sep="")] <- tempprob[,aa] 
                  }
                  
                  ##Kashahara and Shimotsu
                  if (aa != min(alter)) { #to guarantee KS probabilities are within [0,1] bounds
                    pooleddatalistks[,paste("s",wexer,".",aa,sep="")] <- ((tempprobks[,ap])^alpha)*((ml.m.w0ks$probabilities[,ap])^(1-alpha)) #Kasahara and Shimotsu (2012)
                  } else {
                    pooleddatalistks[,paste("s",wexer,".",aa,sep="")] <- 1-apply(((tempprobks[,paste(alter[alter!=min(alter)],sep="")])^alpha)*((ml.m.w0ks$probabilities[,paste(alter[alter!=min(alter)],sep="")])^(1-alpha)),1,sum) #Kasahara and Shimotsu (2012)
                  }
                }
                
                ##Aguirregabiria and Mira (2007)
                if (it==1) {
                  #splitting the previous data into a list 
                  datalist  <- split(pooleddatalist, list(pooleddatalist$PARID))
                  #spatially weighting
                  provweigh <- vector(length=L,mode="list")
                  for(aa in c(1:L)){
                    provweigh[[aa]] <- mapply(weighting,datalist,weigh)
                  }
                  #updating pooled datalist
                  pooleddatalist <- rbindlist(datalist)
                  for(aa in c(1:L)){
                    pooleddatalist[[paste("s",wexer,".",aa,sep="")]] <- unlist(as.vector(provweigh[[aa]]))
                  }
                  rm(provweigh)
                  attr(pooleddatalist,"class") <- "data.frame"
                  
                  #updating expectations
                  expec1 <- subset(pooleddatalist,select=paste("s",wexer,".",alter,sep=""))
                  
                  #Updating pooleddata (long format)
                  for(aa in c(1:L)){
                    pooleddata[pooleddata$alt==aa,][[paste("s",wexer,sep="")]] <- pooleddatalist[[paste("s",wexer,".",aa,sep="")]]
                  }
                  print(paste("AM one step iteration completed",sep=""))
                }
                
                ##Kashahara and Shimotsu
                #splitting the previous data into a list 
                datalistks  <- split(pooleddatalistks, list(pooleddatalistks$PARID))
                
                #spatially weighting
                provweighks <- vector(length=L,mode="list")
                for(aa in c(1:L)){
                 provweighks[[aa]] <- mapply(weighting,datalistks,weigh)
                }
                
                #updating pooled datalist
                pooleddatalistks <- rbindlist(datalistks)
                for(aa in c(1:L)){
                  pooleddatalistks[[paste("s",wexer,".",aa,sep="")]] <- unlist(as.vector(provweighks[[aa]]))
                }
                rm(provweighks)
                attr(pooleddatalistks,"class") <- "data.frame"
                
                #updating expectations
                expec1ks <- subset(pooleddatalistks,select=paste("s",wexer,".",alter,sep=""))
                
                #Updating pooleddata (long format)
                for(aa in c(1:L)){
                  pooleddataks[pooleddataks$alt==aa,][[paste("s",wexer,sep="")]] <- pooleddatalistks[[paste("s",wexer,".",aa,sep="")]] 
                }
                
                if (it>=2) {
                  ndexpectitks <- max(as.matrix(abs(expec1ks-expecks))) # the maximum of all matrix entries #sum(distr(expecks,expec1ks)) #element by element
                }
                print(paste("Max abs diff inner Beliefs KS: ",ndexpectitks,sep=""))
                expecks <- expec1ks
                it <- it+1
              }
              if (it<=maxiter & (ndexpectit < zero | ndexpectitks < zero)) {
                print(paste("Fixed point in beliefs convergence KS reached at iter: ",it-1,sep="")) 
              } else {
                print(paste("Fixed point in beliefs convergence NOT reached at iter: ",it-1,sep=""))
              }
              
              #updating expectations
              expec1o <- subset(pooleddatalist,select=paste("s",wexer,".",alter,sep=""))
              expec1oks <- subset(pooleddatalistks,select=paste("s",wexer,".",alter,sep=""))
              
              if (it0>=2) {
                ndexpectito <- max(as.matrix(abs(expec1o-expeco))) # the maximum of all matrix entries ndJeit <- sum(abs(Jeit[it0,]-Jeit[it0-1,]))
                ndexpectitoks <- max(as.matrix(abs(expec1oks-expecoks))) # the maximum of all matrix entries ndJeitks <- sum(abs(Jeitks[it0,]-Jeitks[it0-1,]))
              }
              print(paste("max Diff Abs ML ",ndexpectito,sep=""))
              print(paste("max Diff Abs ML KS",ndexpectitoks,sep=""))
              expeco <- expec1o
              expecoks <- expec1oks
              
              #To keep the previous probabilities to perform Kasahara and Shimotsu
              ml.m.w0ks <- ml.m.wks
              it0 <- it0+1
            }
            print(paste("Outer ML Heterogeneous iter converged at iter: ",it0-1," using weights: ",wexer,sep=""))
            
            if(wexer == "w") {
              a0ea0[itt0,c(L:(2*(L-1)))] <- summary(ml.m.w)$coeff[paste("(Intercept):",c(2:L),sep="")]
              aea[itt0,c(L:(2*(L-1)))] <- summary(ml.m.w)$coeff[paste("x:",c(2:L),sep="")]
              ded[itt0,c(L:(2*(L-1)))] <-  summary(ml.m.w)$coeff[paste("xw:",c(2:L),sep="")] 
              JeJ[itt0,c((L+1):(2*L))] <- summary(ml.m.w)$coeff[paste("s",wexer,":",c(1:L),sep="")]
              #For KS
              a0ea0ks[itt0,c(L:(2*(L-1)))] <- summary(ml.m.wks)$coeff[paste("(Intercept):",c(2:L),sep="")]
              aeaks[itt0,c(L:(2*(L-1)))] <- summary(ml.m.wks)$coeff[paste("x:",c(2:L),sep="")]
              dedks[itt0,c(L:(2*(L-1)))] <-  summary(ml.m.wks)$coeff[paste("xw:",c(2:L),sep="")] 
              JeJks[itt0,c((L+1):(2*L))] <- summary(ml.m.wks)$coeff[paste("s",wexer,":",c(1:L),sep="")]
            }
            else {
              a0ea0[itt0,c((2*L-1):(3*(L-1)))] <- summary(ml.m.w)$coeff[paste("(Intercept):",c(2:L),sep="")]
              aea[itt0,c((2*L-1):(3*(L-1)))] <- summary(ml.m.w)$coeff[paste("x:",c(2:L),sep="")]
              ded[itt0,c((2*L-1):(3*(L-1)))] <-  summary(ml.m.w)$coeff[paste("xw:",c(2:L),sep="")] 
              JeJ[itt0,c((2*L+1):(3*L))] <- summary(ml.m.w)$coeff[paste("s",wexer,":",c(1:L),sep="")]
              #For KS
              a0ea0ks[itt0,c((2*L-1):(3*(L-1)))] <- summary(ml.m.wks)$coeff[paste("(Intercept):",c(2:L),sep="")]
              aeaks[itt0,c((2*L-1):(3*(L-1)))] <- summary(ml.m.wks)$coeff[paste("x:",c(2:L),sep="")]
              dedks[itt0,c((2*L-1):(3*(L-1)))] <-  summary(ml.m.wks)$coeff[paste("xw:",c(2:L),sep="")] 
              JeJks[itt0,c((2*L+1):(3*L))] <- summary(ml.m.wks)$coeff[paste("s",wexer,":",c(1:L),sep="")]
            }
            if(itt0 == 2) {
              for (aa in c(1:L)){
              corexp <-round(cor(pooleddata[pooleddata$alt==aa,][[paste("s",wexer,sep="")]],pooleddata[pooleddata$alt==aa,]$expw),2)
              pointtx<- (table(pooleddata$alt*pooleddata$choice)[2:(L+1)]/sum(table(pooleddata$alt*pooleddata$choice)[2:(L+1)]))[aa]
              pointty<- median(pooleddata[pooleddata$alt==aa,]$expw)
              cairo_ps(paste(Graphs,"Den3D",Nparish,"Ng",Ng,wexer,aa,".eps",sep=""), width=7, height=7) #to keep transparency from the intervals  
              points(trans3d(pointtx,pointty,max(kde2d(pooleddata[pooleddata$alt==aa,][[paste("s",wexer,sep="")]],pooleddata[pooleddata$alt==aa,]$expw)[[3]]),persp( kde2d(pooleddata[pooleddata$alt==aa,][[paste("s",wexer,sep="")]],pooleddata[pooleddata$alt==aa,]$expw), theta = 30, phi = 30, main=paste("s",wexer,aa," vs ","expw",aa,". corr:",corexp,sep=""), sub = paste("after it:",it0,sep=""), xlab=paste("s",wexer,aa,sep=""), ylab=paste("expw",aa,sep=""), zlab="density",ticktype='detailed')))
              dev.off()
              #For KS
              corexpks <-round(cor(pooleddataks[pooleddataks$alt==aa,][[paste("s",wexer,sep="")]],pooleddataks[pooleddataks$alt==aa,]$expw),2)
              pointtxks<- (table(pooleddataks$alt*pooleddataks$choice)[2:(L+1)]/sum(table(pooleddataks$alt*pooleddataks$choice)[2:(L+1)]))[aa]
              pointtyks<- median(pooleddataks[pooleddataks$alt==aa,]$expw)
              cairo_ps(paste(Graphs,"KS/Den3D",Nparish,"Ng",Ng,wexer,aa,".eps",sep=""), width=7, height=7) #to keep transparency from the intervals  
              points(trans3d(pointtxks,pointtyks,max(kde2d(pooleddataks[pooleddataks$alt==aa,][[paste("s",wexer,sep="")]],pooleddataks[pooleddataks$alt==aa,]$expw)[[3]]),persp( kde2d(pooleddataks[pooleddataks$alt==aa,][[paste("s",wexer,sep="")]],pooleddataks[pooleddataks$alt==aa,]$expw), theta = 30, phi = 30, main=paste("s",wexer,aa," vs ","expw",aa,". corr:",corexp,sep=""), sub = paste("after it:",it0,sep=""), xlab=paste("s",wexer,aa,sep=""), ylab=paste("expw",aa,sep=""), zlab="density",ticktype='detailed')))
              dev.off()
              }
            }
          }
        }
        else {
          print(paste("Outer it: ",itt0," not included due to lack of enough observations in each alternative", sep=""))
          nogooditt0 <- c(nogooditt0,itt0)
        }
      }
      else {
        print(paste("Outer it: ",itt0," not included due to lack of good variation", sep=""))
        nogooditt0 <- c(nogooditt0,itt0)
      }
      
      itt0 <- itt0+1
      nogooditt0 <- unique(nogooditt0)
      gooditt0 <- itt0 - length(nogooditt0)
      print(paste("good worlds:",gooditt0))
    }
    if (length(nogooditt0) >1) nogooditt0 <- c(nogooditt0[2:length(nogooditt0)],itt0:maxitt0) else nogooditt0 <- c(itt0:maxitt0)
    dataJeJ <- as.data.frame(cbind(JeJ,nn0,nnei0)[-nogooditt0,])
    colnames(dataJeJ) <- c(paste("Jexpw",c(1:L),sep=""),paste("Jsw",c(1:L),sep=""),paste("Jswk",c(1:L),sep=""),paste("J",c(1:L),sep=""),"nobs","nnei")
    dataded <- as.data.frame(cbind(ded,nn0,nnei0)[-nogooditt0,])
    colnames(dataded) <- c(paste("dexpw",c(2:L),sep=""),paste("dsw",c(2:L),sep=""),paste("dswk",c(2:L),sep=""),paste("d",c(2:L),sep=""),"nobs","nnei")
    dataaea <- as.data.frame(cbind(aea,nn0,nnei0)[-nogooditt0,])
    colnames(dataaea) <- c(paste("aexpw",c(2:L),sep=""),paste("asw",c(2:L),sep=""),paste("aswk",c(2:L),sep=""),paste("a",c(2:L),sep=""),"nobs","nnei")
    dataa0ea0 <- as.data.frame(cbind(a0ea0,nn0,nnei0)[-nogooditt0,])
    colnames(dataa0ea0) <- c(paste("a0expw",c(2:L),sep=""),paste("a0sw",c(2:L),sep=""),paste("a0swk",c(2:L),sep=""),paste("a0",c(2:L),sep=""),"nobs","nnei")
    
    save(dataJeJ, file=paste(Results,"/dataJeJAll",Nparish,"Ng",Ng,".RData",sep=""))
    save(dataded, file=paste(Results,"/datadedAll",Nparish,"Ng",Ng,".RData",sep=""))
    save(dataaea, file=paste(Results,"/dataaeaAll",Nparish,"Ng",Ng,".RData",sep=""))
    save(dataa0ea0, file=paste(Results,"/dataa0ea0All",Nparish,"Ng",Ng,".RData",sep=""))
    
    #For KS
    dataJeJks <- as.data.frame(cbind(JeJks,nn0ks,nnei0ks)[-nogooditt0,])
    colnames(dataJeJks) <- c(paste("Jexpw",c(1:L),sep=""),paste("Jsw",c(1:L),sep=""),paste("Jswk",c(1:L),sep=""),paste("J",c(1:L),sep=""),"nobs","nnei")
    datadedks <- as.data.frame(cbind(dedks,nn0ks,nnei0ks)[-nogooditt0,])
    colnames(datadedks) <- c(paste("dexpw",c(2:L),sep=""),paste("dsw",c(2:L),sep=""),paste("dswk",c(2:L),sep=""),paste("d",c(2:L),sep=""),"nobs","nnei")
    dataaeaks <- as.data.frame(cbind(aeaks,nn0ks,nnei0ks)[-nogooditt0,])
    colnames(dataaeaks) <- c(paste("aexpw",c(2:L),sep=""),paste("asw",c(2:L),sep=""),paste("aswk",c(2:L),sep=""),paste("a",c(2:L),sep=""),"nobs","nnei")
    dataa0ea0ks <- as.data.frame(cbind(a0ea0ks,nn0ks,nnei0ks)[-nogooditt0,])
    colnames(dataa0ea0ks) <- c(paste("a0expw",c(2:L),sep=""),paste("a0sw",c(2:L),sep=""),paste("a0swk",c(2:L),sep=""),paste("a0",c(2:L),sep=""),"nobs","nnei")
    
    save(dataJeJks, file=paste(Results,"/KS/dataJeJAll",Nparish,"Ng",Ng,".RData",sep=""))
    save(datadedks, file=paste(Results,"/KS/datadedAll",Nparish,"Ng",Ng,".RData",sep=""))
    save(dataaeaks, file=paste(Results,"/KS/dataaeaAll",Nparish,"Ng",Ng,".RData",sep=""))
    save(dataa0ea0ks, file=paste(Results,"/KS/dataa0ea0All",Nparish,"Ng",Ng,".RData",sep=""))
    
  }
}
  sink()
